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Abstract 

We investigate optical Bloch oscillation, Zener tunneling and breathing modes in arrays of optical 
waveguides. We perform a new method of calculation based on the multiple scattering formalism. To 
demonstrate Bloch oscillation and breathing modes, we consider a planar array of parallel cylindrical 
O waveguides with the refractive index gradually varying across the array. We demonstrate that the form of 
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Bloch oscillation may be predicted by means of dispersion law analysis. To demonstrate Zener tunneling, 
we consider a planar array of cylindrical waveguides of two types situated by turn. The band structure of 



■ this array contains two bands separated by a narrow gap. If the refractive indices of waveguides gradually 
O vary across the array, the Zener tunneling leads to the Bloch-Zener oscillation. 



PACS numbers: 
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I. INTRODUCTION 



Nowadays, much attention is devoted to arrays of evanescently coupled optical waveguides which 
are both of the fundamental and practical interest. These arrays are useful in integrated optical 
circuits and other micro- and nanooptical devices, such as optical filters and near-field microscopes. 

The periodic arrays of optical waveguides represent the particular case of low- dimensional pho- 
tonic crystal structures. The general feature of such systems is the existence of photonic band 
structure [1] that is analogous to the electron band structure in solids. Therefore some effects in 
optical lattices may be analogous to some phenomena in solids In this work we consider 

optical counterparts of Bloch oscillation and Zener tunneling. 

Around 1930's it was predicted that an electric field applied to a crystal should induce an 
oscillatory motion of the electrons, known as Bloch oscillation 0, Is]. Besides, in multiband systems 
electrons under an external force can spontaneously transit from one band to another. This effect 
is known as Zener tunneling. 

Optical excitations in arrays of waveguides can perform a similar effects, as it was shown in 
numerous theoretical 6|, |7| and experimental works. The usual pattern to demonstrate the optical 
Bloch oscillation and Zener tunneling is a planar array of parallel waveguides with refractive index 
linearly varying across the array. To produce t he g radual refractive index alteration, one can use 
the thermo-optic SMlO] or electro-optic effects The other pattern is an array of waveguides 
of the same refractive indices^ bu t of different thickness fl2|. Sometimes array of identical gently 



curved waveguides is used 
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The optical excitation coupled into such array propagates along the direction of waveguides 
oscillating in the transverse direction, so the propagation way of the excitation takes the sinusoidal 
form. This phenomena is the optical counterpart of electronic Bloch oscillation in solids. The effect 
of optical Bloch oscillation can be practically used in different optical devices for light steering. 

The optical counterpart of Zener tunneling may take place in presence of two bands separated by 
a gap in the band structure of the array. The superposition of Bloch oscillation and Zener tunneling 



causes the splitting of an optical beam into two breams propagating along different oscillating ways. 



This effect is known as Bloch-Zener oscillation j3, 

In most of works, for theoretical simulation of Bloch oscillation and Zener tunneling the following 
system of equations is used: 
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Here the waveguides are assumed to be directed along the z-axis, j is the number of a waveguide, 
aj{z) is the amphtude of the optical excitation at the j-th waveguide, Pj is the propagation constant 
of the j-th waveguide, 7 is the coupling constant. This system of equations is useful for the 
waveguides of any form, but the parameters 7 and Pj should be obtained experimentally. 

In this paper we use another method of theoretical simulation based on multiple scattering 
formalism (MSF) 16|, |l7|. This method is convenient for the arrays of cylindrical waveguides. Its 
advantage is that the radii and refractive indices of the waveguides are the only data required for 
the calculation, and one has not to obtain any other parameters from an experiment. Besides, 
this method allows to calculate the spatial distribution of electromagnetic field around waveguides 
and inside of them with arbitrary accuracy, as opposed to Eq. ([1]), that allows only to find the 
intensity of optical excitation near every waveguide. 

The MSF is explained in Sect. II. In Sect. Ill we calculate the band structure of a plane array 
of infinite cylindrical rods. The obtained dispersion laws are used in Sect. IV for prediction of 
Bloch oscillation of optical beam in an array of rods with refractive index gradually varying across 
the array. The prediction is confirmed by the direct numerical simulation represented in Sect. 
V. Besides, in Sect. V the so-called breathing mode is investigated. In Sect. VI we investigate 
Bloch-Zener oscillation of optical beam in a plane array of rods of two types situated by turns, with 
gradually varying refractive indices. Finely, in Conclusion we discuss possible practical applications 
of the investigated optical effects and the possibility of further development of method used in this 
paper. 



II. MULTIPLE SCATTERING FORMALISM 



We consider an array of parallel dielectric waveguides directed along the z-axis. We assume 
the waveguides being infinite cylindrical rods. The array is illuminated by a monochromatic wave 
of frequency oj. The velocity of light in free space is supposed to be unit. 

The general idea of multiple scattering formalism is that near the j-th waveguide the incident 
wave can be represented as a linear combination of harmonics with certain values of angular 
momentum m and longitudinal wave vector K: 
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Here r^^^ and (f)^^^ are the polar coordinates of two-dimensional vector r^-') = {x — Xj,y — yj}, and 
the coordinates of the axis of the j-th waveguide. Coefficients pjm{K), qj^{K) are called 
the partial amplitudes of the incident wave. The expressions for functions M^^^^(r), N_^^^^(r) are 
given in Appendix. One can see that (^l^x„^ (r) = 0, therefore harmonics containing functions 
-'^wl'ml^)' •'^ixml^) "^^^ named TM- and TE-harmonics correspondingly. 

The wave scattered by the j-th waveguide can be represented in the similar way, but the 
other functions M^^^^(r), N^^^^(r) enter into the expressions instead of the functions M^^^^(r), 



m=— 00 
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The coefficients ajm{K), hjm{K) are named the partial amplitudes of the scattered wave. The 

('2) i'l) 

functions M;;^„(r), N 

ijjKmi^) ^1^*^ given in Appendix. 
The partial amplitudes of incident and scattered waves satisfy to the following system of equa- 
tions: 



Here k = \J uj'^ — A'^, r^^, are the polar coordinates of two-dimensional vector r^^ = {xj — 
xi, y.j — yi}, Hn{r) is the Hankel function of the first kind, and Sjm{^, K) is the scattering matrix 
for the j-th waveguide. The scattering by a cylindrical waveguide doesn't mix harmonics with 
different longitudinal wave vectors K and with different angular momenta m, but harmonics of 
TEl- and TM-types mix. The formulae to calculate the scattering matrix Sjm{oj, K) are cited in 
Appendix. 



System (jlj) allows to find partial amplitudes of waves scattered by all the waveguides of the 
array. The spatial distribution of field can be calculated by formulae 

N 

E{t, r) = Eincit, r) + Y, ^sca{t, r), 

(5) 

N ^ ' 

H(t, r) = Hi„,(t, r) + ^ H,,,(t, r). 

i=i 

The exact system (jlj) consists of infinite number of equations, containing infinite number of 
variables. The number of equations and variables can be limited, choosing some maximal absolute 
value of angular momentum mmax and taking into account only equations and partial amplitudes 
with m lying in interval — mmax ^ rn < m^^y^. The spatial distribution of field can be calculated 
with required accuracy choosing was demonstrated by the direct numerical 

simulation, that mmax = 2 is enough for qualitative description of optical excitation behaviour. 



III. BAND STRUCTURE CALCULATION. 



Consider the infinite periodic plane array of identical cylindrical waveguides. The array is 
situated in x^-plane, and waveguides are directed along 2;-axis. The distance between two adjacent 
waveguides is a. Below we discuss the eigenmodes of this array and describe the method of band 
structure calculation. 

The system of equations for eigenmodes has the left-hand side coinciding with that in system 
(jl]), and its right-hand side is zero. 

(SU^,K)Y n-^^M - £ 5^e^(«-H^. i/„_(xn,) p^^M = 0. (6) 

Since the array is planar, rij = a\l — j\, (pij = for j > I and (pij = vr for j < I. 
The eigenmodes of an infinite periodical array take the form of Bloch waves characterized by 
the transversal quasi- wave vector k (— vr/a < k < n /a): 

a^m{K) = aUk, K) e*'="^ h,^{K) = b^k, K) (^^^K (7) 
Substituting these expressions to ([6]), a system of equations for a^Sk-, K)^ bm{k, K) is obtained: 
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where 



U^nioJ, k, K) = {Sm{uJ, K)) b^n ' (e"'"' + ("1) 




(9) 



The angular momentum takes values — m^ax < m < mmaxj so ([8]) is a homogeneous linear 
system of 4mmax + 2 equations with the same number of variables. If this system is represented in 
matrix form, its matrix U{u, k, K) is composed of (2minax + 1)^ matrices Umn{^, k, K). The system 
([8]) possesses a nontrivial solution when the matrix U{uj, k, K) is singular: det U{(jJ, k, K) = 0. 

Using the technic described above, we calculated the band structure for an array of cylindrical 
rods of unit radii {R = 1) made of GaAs (refractive index Ur = 3.5). The rods are situated next to 
each other, so the period of the array is a = 2R = 2. We calculated the dependence of longitudinal 
wave vector K on transverse quasi-wave vector k for a fixed frequency u = 0.77r/a (below we will 
use the term "dispersion law" for the dependence K{k)). The approximation mmax = 2 was used. 
It was shown by the direct numerical simulation, that for the chosen m^ax the dispersion law K{k) 
can be found accurate within 2%. The dispersion curves K{k) are presented at Fig. [H 

In Fig. [T] several dispersion curves corresponding to several different bands are illustrated. 
Below we consider one of the bands, that is noted by letter "A". This band is convenient for further 
investigation, since it doesn't overlap with other bands. 

Since the parameters K and k of eigenmodes are connected by the dispersion laws, one of 
arguments in notations am{k,K), bm{k,K) for partial amplitudes is unnecessary, so below the 
partial amplitudes are denoted am.{K), hm{K). 

IV. BLOCK OSCILLATION PREDICTION ON BASIS OF DISPERSION LAW. 

Below we consider the Gaussian beam propagating in the array of waveguides. The partial 
amplitudes describing this excitation are represented by formulae 



-r2 {K-Kof 
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Figure 1: Dispersion curves K(k). 
In this case the field distribution takes the form 



E{t, r) = e-*^*u(r) exp I - ^^^^ +ikoX + iKoz 



At 

H(t, r) = e-'"V(r) exp <| -i^/^-^ + ikoX + iKoz 



Here u(r), v(r) are the functions periodically depending on x, ko is connected with Kq by the 
dispersion law, Kq = K{kQ), and v = dK/dk{ko). The formulae ( ITTj) are correct for enough large 
values of r. 

It follows from Eqs ( ITTj) . that in the periodical array of identical waveguides the optical excita- 
tion propagates along the straight line x{z) = vz, and the direction of propagation is defined by 
the dispersion law K{k). 

But the situation changes dramatically, if the optical characteristics of waveguides (such as 
thickness or refractive index) gradually vary across the array. 

One can mentally divide the array to sections much wider than the optical beam, but enough 
narrow for one could assume the waveguides into a section to be identical. One can attribute 



a local dispersion law K{k) to every section. Therefore, the direction of the beam propagation 
should be different in different sections, and the propagation way of optical excitation should be 
curved. A certain form of the propagation way can be predicted by calculating the dispersion law 
for arrays with different refractive indices of the waveguides. 

For example, consider an array of = 100 cylindrical rods of unite radii. The refractive index 
of a rod in the middle of the array (j = 50) is nj!° = 3.5, and the difference between the refractive 
indices of two adjacent waveguides is nl — nl'^^ = 0.01, i. e. 

ni = 3.5 - 0.01 (j - 50) (12) 

The section at the middle of the array is similar to the array considered in the previous section. 
So, we choose the parameters of optical beam according to the dispersion curve represented in 
Fig. [Hand marked by letter "A". The band corresponding to that dispersion curve lies in the range 
2.35 < K < 2.78, so we choose Kq = 2.565 exactly at the middle of the band. The frequency of 
the excitation u = 0.357r. 
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Figure 2: Dispersion laws for different refractive indices. The straight line corresponds to K = 2.565. 



If the refractive index changes, the dispersion curve "A" shifts, as it is shown in Fig. [2l We have 
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found that the longitudinal wave-vector Kq lies into the band "A" if the refractive index varies in 
the interval 3.35 < rir < 3.65. Therefore, the optical excitation can propagate in a part of the 
array where the refractive indices of waveguides belong to the mentioned interval, i.e. between the 
35-th and 65-th waveguides. 

So, we have predicted the amplitude of Bloch oscillation. But one can also predict the period of 
Bloch oscillation and the way of optical beam propagation. For this purpose, the value v = dK/dk 
for Kq for different values of refractive index rir should be calculated. For the obtained dependence 
the notation f (n,.) will be used. The way x{z) of optical beam propagation is determined by the 
differential equation 




(13) 



where the function nr{x) is obtained by the interpolation of dependence of the waveguide 
refractive index nl on the waveguide number j: 

ririx) = 3.5 - 0.01 (^^ - 50) . (14) 

Eq. (fT3|) can be integrated numerically. The way of optical beam propagation, obtained from 
this equation, has the form of periodical oscillation, as represented in Fig. [31 The period of the 
obtained oscillation is Az ~ 220a = 440 (remind that a = 2). 



V. DIRECT CALCULATION OF BLOCH OSCILLATION AND BREATHING MODE. 

In this section we represent the results of direct calculation of Gaussian beam propagation, 
based on numerical solution of Eq. |H We consider the same array as in the previous section. The 
array is illuminated by an incident wave that is defined by partial amplitudes 

pUK) = Pm exp L^^i^^M. + ,koa {j - jo) | , 

qjUK) = e-^' (^-^0)^ exp ^\J'^ + ik,a {j - j,) j . 

This incident wave illuminates the finite area of the array. The parameters a and r define the 
width of the illuminated area along x-axis and z-axis correspondingly. We take jo = 50, i. e. the 
incident wave illuminates the middle of the array. The parameters Kq and are connected by 
the dispersion law marked by letter "A" in Fig. [H We take Kq = 2.565, exactly at the middle of 
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Figure 3: Optical excitation propagation way obtained by the analysis of dispersion law. 

the band "A" (2.35 < K < 2.78), and the corresponding = 0.989. The parameter r is chosen 
so that the peak of the function e"'^^^^"^")^ fits into the band "A": r = 6/AK ^ 14, where AK 
is the width of the band "A". The parameter a = f (/cq) t, where f (/cq) is defined by the dispersion 
law: f (fco) = dK/dk{ko). Here v{ko) ^ 0.5, so cr = 7. 

The parameters p^, Qm are chosen so that the incident wave excites the eigenmodes of the 
array effectively. For our calculation we chose go = 1 and all the other pm, Qm are zeros. The 
computation is performed for the approximation mmax = 2. 

The result of the computation is presented in Fig. HI As expected, the obtained way of optical 
beam propagation has a periodical form. It oscillates between the 35-th and 65-th waveguides, 
and the period of oscillation is Az = 440. It is remarkable that the form of oscillation obtained 
by the numerical solution of Eq. (jll) coincides exactly with that obtained by the dispersion law 
analysis. 

Besides the Bloch oscillation, we consider the so-called breathing mode jo, Q, [l^. Such kind of 
optical excitation arises when only one waveguide of the array is illuminated by the incident wave. 



10 



3000 



2500 



2000 

■M 1500 
O 

u 

1000 



500 




10 20 30 40 50 60 70 80 90 100 
Number of a waveguide j 



Figure 4: Optical excitation propagation. 

The characteristic feature of breathing mode is the periodical spreading and focusing behaviour. 

We assume that the incident wave illuminates a short section around z = of the 50-th 
waveguide situated at the middle of the array. To simulate this situation, we take the partial 
amplitudes of the incident wave as follows: pjm{K) = 0, qjm{K) = 0, qjo{K) = 1. The result of 

m^O 

direct numerical calculation is represented in Fig. O 

VI. BLOCH-ZENER OSCILLATION. 

The optical Bloch-Zener oscillation in an array of optical waveguides can take place if the 
band structure consists of several bands separated by gaps. If the refractive index of waveguides 
gradually varies across the array, the band structures of two different sections of the array are 
shifted relative to each other. Therefore, the lower band of one section can overlap the upper band 
of another section. So, the optical beam can partially tunnel from one section to another. This 
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Figure 5: Breathing mode. 

leads to that the optical beam divides into two beams propagating along two different oscillating 
ways. 

To demonstrate this effect, it is convenient to consider an array with the band structure contain- 
ing two bands separated with a narrow gap. We take the array of waveguides of two types situated 
by turns. The refraction indices of waveguides of the first and the second type are riri — 3.5 and 
nr2 — 3.55 respectively. The band structure of this array contains two bands 2.38 < K < 2.48 and 
2.54 < K < 2.8 (we suppose the frequency cu = 0.357r, as in previous sections). 

We introduce a small variation of refraction indices of waveguides: 

= 3.5 - 0.005 (j - 50) for odd j, 

(16) 

= 3.55 - 0.005 (j - 50) for even j. 
The array is illuminated by the incident wave defined by the partial amplitudes pjm{K), qj^{K), 



12 



that are given by formulae (fT5|) . The parameters Kq, ko, r, a entering to these formulae are 
chosen according to the principle similar to that described in the previous section. The parameter 
Kq = 2.67 is taken exactly at the middle of the band 2.54 < K < 2.8, the parameter = 0.44 is 
connected to Kq by the dispersion law K{k), the value of r = 23 is taken so that the peak of the 
function e"^"^ (^"-^o)^ entirely fits into the band 2.54 < K < 2.8, and a = dK/dk{ko) r = 10. 
The result of the calculation is represented on Fig. [6l 
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Figure 6: Bloch-Zener oscillation. 

VII. CONCLUSION. 

In this paper three phenomena are considered — Bloch oscillation, Bloch-Zener oscillation and 
breathing modes in planar arrays of optical waveguides with gradually varying refractive index. 
We suggest a new method to investigate this subject, based on the multiple scattering formalism. 
This method has several advantages over the traditional method based on Eq. ([1]). The MSF 
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allows to find the spatial distribution of field with any required accuracy, while the traditional 
method gives only the intensity of optical excitation. Besides, the input data for MSF are the 
geometrical properties of the array and refractive indices of waveguides, while the traditional 
method requires some data that should be obtained experimentally, such as the longitudinal wave 
vectors of eigenmodes of waveguides and coupling constants. 

The MSF represented in this paper is convenient only for the waveguides of cylindrical form, 
because in this case the scattering matrix can be calculated easily. However, this method can 
be applied for the waveguides of another shape, but in this case it would be more difficult to 
calculate the scattering matrix. Besides, the scattering by noncylindrical waveguides would mix 
the harmonics with different angular momenta. So, if the shape of the waveguides is enough 
complicated, one should take into account the harmonics with enough high angular momenta, and 
the calculation would be difficult. At the same time, for the cylindrical waveguides it is enough to 
take into account the harmonics with |m| < 2, as it is shown in this work. 

The considered phenomena may be useful for different optical applications, such as steering, 
splitting, focusing and defocusing of hght. The method represented in this work allows to produce 
the numerical simulation without need of experimental investigation of components of optical 
devices. 
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APPENDIX. 

1. Formulae for functions M[j'^^(r) and N[^^^^(r): 




iK 



(17) 




(18) 




Here x = y/co'^ — K"^, Jm{x) is Bessel function. 



2. Formulae for functions M, 




(r) and N| 




(r): 




iK 



(19) 
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Here Hm{x) is Hankel function of the first kind. 
3. Formulae for scattering matrix: 

Consider an infinite dielectric rod situated along the 2;-axis. The radius of the rod is R, and 
its refractive index rir- It is illuminated by a monochromatic wave of frequency ou with certain 
longitudinal wave vector K and angular momentum m. This wave is defined by two partial 
amplitudes pm{K), q^^K). The scattered wave possesses the same frequency cj, longitudinal wave 
vector K and angular momentum m. It is defined by partial amplitudes am{K), hm{K). 

Partial amplitudes of incident and scattered waves are connected by the scattering matrix 
Sm{uJ,K): 

hm{K) ) \qm{K) ) 

To formulate the expression for matrix Sm{i^, K), we introduced some notations: 

x^Vco^-K^, e = a^K/2x, ^ = a;/2x, ^^^^ 



U- ^ Hm-l{0 - Hm+liO^ Mo = -f^m(O) U+ = Hm-l{0 + ^m+l{0 ^ (23) 

Mn = I , M21 = ' 



au^ iPu-l \ —uq 



M„= I I. M.= r" - 1, (24) 

-ai v+ -i(3i V. 

iaw- —Bw+ 
—awj^ —i/3w^ 
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Using the introduced notations, we write down the expression for matrix Sm{i^, K) 
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